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ABSTRACT 

We have developed ID time-dependent numerical models of accretion discs, using an 
adaptive grid technique and an implicit numerical scheme, in which the disc size is 
allowed to vary with time. The code fully resolves the cooling and heating fronts 
propagating in the disc. We show that models in which the radius of the outer edge 
of the disc is fixed produce incorrect results, from which probably incorrect conclu- 
sions about the viscosity law have been inferred. In particular we show that outside-in 
outbursts are possible when a standard bimodal behaviour of the Shakura-Sunyaev 
viscosity parameter a is used. We also discuss to what extent insufficient grid resolu- 
tions have limited the predictive power of previous models. We find that the global 
properties (magnitudes, etc. ...) of transient discs can be addressed by codes using a 
high, but reasonable, number of fixed grid points. However, the study of the detailed 
physical properties of the transition fronts generally requires resolutions which are out 
of reach of fixed grid codes. It appears that most time-dependent models of accretion 
discs published in the literature have been limited by resolution effects, improper outer 
boundary conditions, or both. 
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1 INTRODUCTION 

The thermal- viscous accretion-disc i nstabil ity model is more 
than 15 years old (see Cannizzo (1993b) for a historical 
overview). It is widely accepted that it provides the cor- 
rect description of dwarf-nova outbursts and of ('soft') X-ray 
transient events. When, however, observations of these sys- 
tems are compared with predictions of the model, the a gree- 
ment is far from perfect (e.g. Lasota fc Hameury (1998) and 



mass-transfer mod el were thought to contradict observations 
( |Cannizzo, 1993b ) and the outer disc radius behaviour dur- 
ing and after outbursts seemed t o favour the disc instability 
model (tchikawa & Osaki, 1992). The demise of the mass- 



transfer instability model was, however, caused by the lack 
of a physical mechanism which would trigger it. 

With one model left it became important to establish 
just what its predictions are, and not merely whether it is 



better (or worse) th an the competing model (Pringle, Ver- 



references therein). It is sometimes also unclear what the bunt & Wade 1986). The first systematic stud y of th e disc 



predictions of the model are. One of the reasons for these 
uncertainties is the existence of various, different, versions of 
the model. From the very beginning these versions of the disc 
instability model differed in assumptions about viscosity and 
boundary conditions; they differed in the amount of mat- 
ter accreted during the ou tburst, the shapes of light-curves, 
etc. (see Cannizzo (19931)). At that time these differences 



instability model was presented by Cannizzo (1993a), who 
analysed the importance of various terms in the disc evolu- 
tion equations and the influence of the numerical grid resolu- 
tion on t he out burst properties. Ludwig, Meyer-Hofmeister 
& Ritter ( 199^ ) studied general properties of disc outbursts, 
such as the location of the insta bility that triggers them. Re- 
cently Ludwig & Meyer (1998) analysed non-Keplerian ef- 



seemed to be less important than the differences between 
the disc instabilit y model and the comp eting, mass-transfer 
instability model (Bath & Pringle, 1981). The exponentially 
decaying tails of theoretical light curves predicted by the 



fects which may arise during front propagation. The general 
conclusions of this group of studies were that non-Keplerian 
effects are negligible, that a few hundred grid points pro- 
vide a sufficient resolution for the calculation results to be 
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independent of the number of grid points, and, finally, that 
with the usual assumption (see Smak (1984b) of a jump in 
the value of the viscosity parameter a, the model produces 
only 'inside-out' outbursts, i.e. outbursts starting in the in- 
ner disc regions. This last conclusion, if true, would entail 
changing the standard viscosity law (in which the a param- 
eter is constant in the hot and cool branch of the E - T e g 
curve) because outburts starting in the outer disc regions 
are clearly obse rved in classical dwarf-nova system SS Cyg 
( Mauche, 1996). This law a lready had to be modified when 
it was found ( Smak, 1984b ) that in order to get lightcurves 
similar to those observed in dwarf novae, a in outburst had 
to be larger than a in quiescence. 

The absence of 'outside-in' outbursts in these studies, 
however, is just the result of keeping the outer disc radius 



Os- 



constant i n the calculat ions (Section 4.1) (Ichikawa 

pmak, 1984b[); from this point of view, there is 



aki 

no reason to modify the viscosity prescription. This does 
not in itself prove that changes in viscosity are correctly de- 
scribed by the bimodal behaviour of the a-parameter (see 
e.g. Gammie & Menou 1998), but the reasons given for pre- 
ferring other versions (the exponential decay from outburst 
being the principal one) are not compelling, and these ver- 
sions involve more fundamental changes in the disc physics 
(see Lasota & Hameury 1998 for discussion a nd ref erences). 

For example, Cannizzo, Chen & Livio ( 1995 ) use the 
formula a = ao (H/R) n , but to make the model work they 
have to 'switch off' convection. It is interesting, therefore, to 
recall that Faulkner, Lin & Papaloizou (1983) found dwarf 
nova outburst with a c onstant, but their model was criti- 
cized dCannizzo, 1993b| ) because they claimed that convec- 
tion has only a minor influence on the energy transport in 
the disc. These are not formal problems because the constant 
a models predict optically thin quiescent discs, whereas in 
bimodal a models the quiescent disc is optically thick. There 
is observational evidence that dwarf nova discs in quiescence 
are optically thin (see Home, 1993 and references therein). 

Conclusions about the number of grid points required 
to get resolution-invariant results seem, on inspection, too 
optimistic, especially because fronts are not resolved, a point 
which is particularity worrying for the heating fronts. 

The present situation of the disc instability model seems 
to be confused. Various versions are based on different as- 
sumptions about the physical processes in the disc and nu- 
merical codes suffer either from incorrect boundary condi- 
tions or from insufficient resolution or from both. Quite of- 
ten, in the case of explicit codes the resolution is limited by 
the required computer time. 

In this article we describe a numerical model of time- 
dependent accretion discs, using an adaptive grid technique 
and an implicit numerical scheme, in which the disc size 
is allowed to vary with time. This numerical scheme allows 
rapid calculations of disc outburst cycles at very high resolu- 
tion. These properties allow an easy comparison with other 
versions of the model and a systematic study of its various 
assumptions. 

In the near future we will use our code to model var- 
ious properties of dwarf novae and X-ray transients. The 
model was ahead used to model properties and outbursts 
of the dwarf-nova WZ Sge (|Lasota, Hameury fc Hure, 1995 



Hameury, Lasota & Hure, 1997) an d the rise to outburst o f 
the X-ray transient GRO J1655-40 ([Hameury et al., 1997 ). 



In §2 we discuss the time-dependent equations describ- 
ing the disc radial structure and the implicit method used 
to solve them with a high spatial resolution. The vertical 
structure of the disc, and hence the heating and cooling 
terms that enter the time-dependent energy equation, are 
considered in §3. In §4 we present the results of our calcu- 
lations and we discuss the importance of having sufficient 
numerical resolution and a correct boundary condition at 
the outer edge of the disc. 



2 TIME-DEPENDENT ACCRETION DISCS 
2.1 Disc equations 

The basic equations for mass and angular momentum con- 
servation in a geometrically thin accretion disc can be writ- 
ten as: 
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where E is the surface column density, M ex t(r) is the rate 
at which mass is incorporated into the disc at point r, v T 
the radial velocity in the disc, j = (GMir) 1 ^ 2 is the spe- 
cific angular momentum of material at radius r in the disc, 
Q K = (GMi/r 3 ) 1/2 is the Keplerian angular velocity (Mi 
being the mass of the accreting object), v is the kinematic 
viscosity coefficient, and j'k the specific angular momentum 
of the material transfered from the secondary. Tt\& is the 
t orque due to tidal forces, for which we use the value of Smak 
(19841), derived from the d etermination of tidal torques by 
Papaloizou & Pringle (1977)): 



T t id = ctorvY, 



(3) 



where u) is the angular velocity of the binary orbital motion, 
c is a numerical coefficient taken so as to give a stationary 
(or time averaged) disc radius equal to some chosen value, 
and a is the binary orbital separation. 

Equation ^ is the standard form of angular momen- 
tum conservation, in which it is assumed that the azimuthal 
veloci ty v$ has the Keplerian value Vk- Ludwig & Meyer 
have considered deviations from Keplerian motion 



([l99§) 



arising from strong pressure gradients (i.e. in the heating 
and cooling fronts) and found that these deviations only 
have a marginal influence on the outburst behaviour. 
The energy conservation equation is taken as: 



dt 



2(Q+ -Q- + J) 5RT C 1 d{rv, 



CpE 



jiCp r dr 



where Q + and Q~ are the surface heating and cooling rates 
respectively. They are usually taken as Q + = (9/8)i/Efijc 
and Q~ = <jT^ a , T e s being the effective temperature (e.g. 
Cannizzo ( 1993a )). As described in §3 we use slightly differ- 
ent forms of these rates. The term J accounts for the radial 
energy flux carried either by viscous processes or by radia- 
tion. In the following, we neglect the radiative flux which, 
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according to Ludwig fc Mey er ( 199£ ) , is negligible, while ac- 
cording to Cannizzo (1993a) it is comparable to the viscous 
radial flux. 

If viscosity is due to turbulence, J can be estimated in 
the framework of the a parametrization. The flux carried in 
eddies with characteristic velocity v e and size l e , is: 



3:1 re sonance (see e.g Frank, King & Raine (1992)). Osaki 



or 

in which case 

J = l/rd/dr(rF c ) 



3 0T C 



(5) 



(6) 



A similar expression is obtained for radiative transport. Note 
that t his exp ression differs slightly from that used by Can- 
nizzo (1993a) which cannot be written as the divergence of 
a flux. Other prescriptions for J exist, in time-dependent 
simulations- they give results very similar to those obtained 
using Eq. (|). 



2.2 Boundary conditions 

So far there is no clear understanding of the physics of the 
interaction between the disc and the stream (see e.g. Ar- 
mitage & Livio, 1997), and a fraction of the stream may 
spill over the disc. A proper treatment of the precise way 
in which matter is incorporated into the disc is far beyond 
our scope, and we use here the simplest, but very reasonable 
assumption that mass addition at the outer edge of the disc 
occurs in a very narrow region, so that the disc edge is very 
sharply defined. We can then write M ex t(r) = MtrS(ro(t)— r) 
and E = T.oY(ro(t) — r), where Mtt is the mass transfer rate 
from the secondary star, Y is the Heavyside function, S is 
the Dirac function and Eq, the surface column density, is a 
smoothly varying function. The cancelation of the 5 terms in 
the equation for mass and angular momentum conservation 
yields two boundary conditions, which can be written in the 
form: 



M tr = 27rr£o(ro — «r,o 
and 

, ! - 



1 



37T!/Eo, 



(7) 



(8) 



where the index denotes quantities measured at the outer 
edge, and rk is the circularization radius, i.e. the radius at 
which the Keplerian angular momentum is that of the mat- 
ter lost by the secondary star, and ro is the time derivative 
of the outer disc radius. It is worth noting that the presence 
of a torque Tya is required in this formulation, and it can 
be easily seen that no steady solutions exist when Ttid = 0. 

Conditions (Q, ^) take into account the fact that the 
outer edge of the disc can vary with time, and its position 
is controlled by the the tidal torque Ttid- Variations of the 
disc radius are obs e rved during outburs t cycles in dwarf no- 
(|Smak 



1984a 



O'Donoghue, 1986 



Wood et al., 1989 



Wolf et al., 1993), with a rapid rise during the outburst 



and a decrease of the order of 20 percent during decline; 
such variations have been considered as a strong argument 
in favour of the disc instability model for dwarf novae, as op- 
posed to the mass transfer instability model. The presence of 
'superhumps' during SU UMa's superoutbursts is explained 
by the disc radius becoming larger than the radius for the 



( [1996D assumes that the superoutbursts themselves are due 



to a 'tidal-thermal instability' which sets in when the disc's 
outer radius reaches the 3:1 resonance zone; see however 



Smak ( |1996| ). 

Nevertheless it is often assumed, when modeling the 
normal outbursts, that the outer edge of the disc is fixed at 
a given radius, in which case Eq. ([j]) is used with ro = 0, 
and Eq. (^) is replaced by r = ro; the tidal torque Ttid is 
also neglected in Eq. (Q). This is equivalent to assuming that 
Ttid is negligible at r < ro and becomes infinite at r = ro. 
In view of the steep functional dependence of Ttjd(r), this 
might seem a reasonable approximation; however, we shall 
show that this is not the case at all, and that results obtained 
with both boundary conditions differ very significantly. An 
interm ediate formulation has been proposed by Mineshige & 
Osaki (1985), who assume that the viscous stresses vanish at 
a fixed outer radius. This enables matter carrying angular 
momentum to leave the disc at its outer edge, at a rate which 
is comparable to the mass transfer rate. 

We take, as usual, E = at the inner edge of the disc 
(more exactly, E equal to a small value). 

The thermal equation being a second order partial dif- 
ferential equation in r, two boundary conditions are re- 
quired. However, except across a transition front between a 
hot and a cool region, the dominant terms in Eq. (Q) are Q + 
and Q~ . The highest order terms in the thermal equation are 
therefore negligible in almost all of the disc; the solutions of 
such equations are known to develop boundary layers, which 
adapt the internal solution (in which the highest order terms 
are neglected) to the boundary conditions. These boundary 
conditions are thus of no physical importance, and should 
be such as to minimize numerical difficulties. We have here 
taken dT c /dr = at both edges of the disc. 



2.3 Numerical method 

We solve the set of equat ions (|l| [], ^) using a method de- 
scribed by Eggleton ( 1971 ). This method, which uses a vari- 
able mesh size, is a g eneralization of the Henyey method 
(Henyey et al., 1959), designed to solve the set of non- 



linear equations describing the internal structure of stars. 
The solutions of these equations have steep gradients, both 
at their surface and in the vicinity of the thin burning shells 
which appear after their evolution off the main sequence; 
moreover, the position of the shells is not fixed, but may 
vary relatively rapidly across the stellar envelope. 

The natural abscissa, r in our case, is considered as an 
unknown variable of a new parameter q in the range 0-1 
over the grid. The variable r varies according to: 



dr 
dq 



$x W(r,T c ,£, 



(9) 



where W is a function which becomes small when the 
radial derivatives of T c or E become large, and $ is a nor- 
malization constant which adjusts the grid to the physical 
range covered by the variable r, so that 



Here, we take 



(10) 
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1 + 0.05 



d\n(uY,r 



l/2\ 



d In r 



0.1 



<91nT c 



91nr 



(11) 



The thermal and viscous equations ([jj ^) are equivalent to 
a set of four first order differential equations. These equa- 
tions plus the two equations [m]) defining the grid are 
discretized and solved using a generalized Newton method. 
They can be written in the form: 



dq 



Gi(f), (< = !,/), 



(12) 



where I = 6 is the number of equations, the f J are the 6 
variables r, $, T c , E, d(vY,r 1/2 )/dr and dT c /dr, and the d 
terms contain the time-derivatives. The discretized equation 
takes the form: 

Hfl) - Fi(fLi) ~ HPiG(fi) + (1 - fr)G(fU)] = 0, (13) 

where 5q = 1/(N — 1) is the mesh size, f 3 k is the value of 
the quantity f J at the grid point k (k = 1,N), and the 
Pi are weights in the range [0, 1]. Values of Pi equal to 0.5 
lead to a second order scheme. The set of equations ( |l3| ) is 
linearized (including Eq. (pf|ic|) by numerical differentiation. 
This method has the advantage of being fully implicit, but 
convergence may be a problem when the initial guess is not 
close enough to the actual solut ion. 
Whereas in Eggleton ( 1 



1971 



) the time derivatives are ex- 
pressed in their Eulerian form, we directly estimate the La- 
grangean values by using either a linear or a cubic spline 
interpolation of the quantities calculated at the previous 
timestep, on the previous grid. The latter is more accurate, 
but requires a higher number of points when a transition 
front is present in the disc in order to avoid numerical oscil- 
lations. 

This formulation does not guarantee that the disc mass 
is a conserved quantity within roundoff errors. Instead, the 
disc mass has to be explicitly calculated by integrating E 
over the disc extension, and deviations from mass conser- 
vation (easy to calculate given the inner and outer mass 
accretion rates in the disc, together with the time step of 
integration) provide valuable information on the quality of 
the code. We checked that all adaptive grid models presented 
here do indeed conserve mass accurately, except when we use 
too small a number of grid points (typically 100 or less). In 
the latter case, convergence is often a problem. 

In order to improve the stability of the numerical 
scheme, we have used values of Pi equal to for the two 
equations defining the derivatives of T c and E, and to 1 for 
the two equations defining their second derivatives when the 
radius of the outer edge of the disc is kept fixed; we have 
also calculated the Lagrangean derivatives using linear in- 
terpolations. However, when ro is allowed to vary, mass con- 
servation is more difficult to enforce; we have taken all Pi 
equal to 0.5 in that case, and used cubic splines to calculate 
the time derivatives. The timestep is chosen on empirical 
grounds and depends on the quality of the previous conver- 
gence. We find that initiating the disc with a hot steady 
solution or a globally cold state with an arbitrary E profile 
gives, after a short transient period, equivalent results. 



3 DETERMINATION OF THE EFFECTIVE 
TEMPERATURE 

The vertical structure equations are very similar to those 
describing the internal structure of a star, with the notable 
exception that energy is dissipated everywhere in the verti- 
cal structure, including the optically thin regions, and that 
the vertical gravity comes from the central object. 

In order to follow the time evolution of accretion discs, 
one needs to know the thermal imbalance (Q + — Q~) as a 
function of the central temperature T c and the surface col- 
umn density E at any radius in the disc. One should then 
be able to solve the full 2D energy transfer problem, which, 
in the thin disc approximation, is assumed to decouple into 
radial (Eq. (^j)) and vertical equations. The latter, together 
with the reasonable ass umption of ve rtical hydrostatic equi- 
librium, are as follows (3mak, 1984b): 



dP „2 



dz 
dlnT 
dlnP 



= 2p, 



= V, 



dz 2 dz 



(14) 
(15) 
(16) 
(17) 



where P, p and T are the pressure, density and temperature 
respectively, ? is the surface column density between —z 
and +z, g z = Q^z the vertical component of gravity, F z the 
vertical energy flux and V the temperature gradient of the 
structure. This is generally radiative, with V = V ra d, given 
by: 

kPFz 

V rad - -T= , (18) 

4^radCg z 

P r ad being the radiative pressure. When the radiative gra- 
dient is superadiabatic, V is convective (V = V CO nv). The 
convective gradient is calculated in the mixing length ap- 
proximation, with a mixing length taken as H m \ = a m iHp, 
where Hp is the pressure scale height: 

HP= P9* + (PpW (19) 

which ensures that Hp is smaller than the vertical scale 
height of the disc. We have (Paczyhski, 1969): 



Vad + (V rad - V a d)V(V + A) 



(20) 



where V a d is the adiabatic gradient, and Y is the solution 
of the cubic equation: 



9 



J^—Y 3 + VY 2 + V 2 Y- 



V = 



(21) 



where r m i = npH m i is the optical depth of the convective 
eddies. The coefficient V is given by: 



3 + g 2 z HjyCp / dlnp - 

3r m i ) 512a 2 T a H P \d\nT / P 



( dlnp \ 
KdlnTJi 

X(Vrad - Vad) 



(22) 



Here we take a m i = 1.5, which is appropriate for solar type 
stars (ami ra nges fr om 1 to 2 in solar models, see e.g . Guz ik 
& Lebreton (1991) and Demarque and Guenther (1991)). 
One should note that a m ; = 1 is generally used in the lit- 
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crature; as t he disc vertical s tructure depends on the as- 
sumed ami (Cannizzo, 1993b), we expect some difference 



-aQ.K 



from previous results. Ft is a time-dependent contribution 
that includes terms resulting from heating/cooling and con- 
traction/expansion. This contribution is basically unknown, 
but it is most often assumed as proportional to P; if this is 
the case, then Eq. dl7h can be replaced by: 



dz 



-oi c sQ.kP, 



(23) 



where a c g may be considered as some effective viscosity. 
This coefficient is not known a priori and is generally differ- 
ent from the true viscosity parameter; the difference a — a e g 
measures the departure from thermal equilibrium. The ver- 
tical structure, and hence the disc effective temperature at 
any given point, can therefore be obtained assuming that 
the disc is in thermal equilibrium, but with some unspeci- 
fied a e ff different from a. 

Even though this approximation is not unreasonable, 
and can be shown to be valid in the simplest case, th at of a 



perfect gas and homologous cooling and contraction (Smak 
1984b|j it is by no means justified and constitutes a signif- 
icant drawback to the model. Other approximations have 
been used in the past; it is for example possible to assume 
that F z has a given vertic al profile; this is the approach of 
Mineshige & Osaki (1983) who assumed that F z increases 
linearly with altitude up to a saturation value. This is not 
drastically different from assuming that F t is proportional 
to pressure, since in both cases, energy is released in the 
densest regions of the disc, i.e. within one or a few scale 
heights from the mid plane. We have compared the results 
given in both cases, and we found that, for a given (E,T C ), 
T c being the central temperature, the effective temperatures 
obtained agreed to within 5 to 10 %. Therefore, the cooling 
fluxes in time-dependent accretion discs are not determined 
better than to within 50%. 

It would be highly desirable to have a better determi- 
nation of Q~; however, this would require at least the in- 
troduction of explicitly time-dependent terms in the vertical 
structure equations, and that would couple the radial and 
vertical structure of the disc. This would enormously in- 
crease the required amount of computing time, making it 
very difficult to follow a complete outburst cycle with a rea- 
sonable number of radial grid points. 

The viscous dissipation term 3/2aQ,i(P in Eq. ( |l7| ) is 
a direct consequence of the assumption that the radial- 
azimuthal compone nt of the viscous stress t ensor can be 
written T rt j, = —otP (Shakura & Sunyaev, 1973). One should 



also note that the a prescription is used here in its local 
version, i.e. the local energy dissipation is assumed to be 
proportional to the pressure, which is by no means obvious. 

Our approach is to determine Q~ for a number of values 
of T c , E, and r and then perform a linear interpolation. 
In what follows, we have taken 70 equally logarithmically 
spaced values for r (190, 120 for T c , E respectively). This 
may seem time-consuming, but once the initial grid has been 
calculated, it is much more precise than the approach in 
which a sparse grid is used to construct analytic fits of Q~ . 
The E — Teg "S-curves" are found numerically by looking 
in the grid for values of E, T c , a and r which satisfy the 
thermal equilibrium condition. 

In thermal equilibrium the heating term Q + is given by: 



Pdz, 



(24) 



which is close, but not exactly equal, to the heating term 
(9/8) ^Ef2^ generally used in the literature. Setting v — 
2/Za.cl /Qk, and assuming that c s is constant, one recovers 
the classical expression. This slight difference might appear 
as a minor inconsistency, but it may lead to unwanted nu- 
merical effects, since the equilibrium situation found when 
solving the radial structure would be different from that 
found in this section, and the disc may become unstable for 
a surface density slightly different from the expected critical 
value. For this reason, we prefer to use the exact value of the 
integral in (|2^) rather than its approximation; it is stored 
and interpolated in order to provide an accurate Q + to the 
time- dependent runs. 

We solved the vertical structure equations in two dif- 
ferent approximations for the radiative transfer of energy, 
namely in the optically thick approximation and in the grey 
atmosphere approximation. 



3.1 Optically thick case 



The set of equations ( |14| - fl6| , |23j) is integrated between the 
disc midplane and the photosphere, defined as the point 
where: 

2 



krP 



(25) 



where is the Rosseland mean opacity. 

The other boundary conditions are z = 0, F z — 0, T = 
T c , ? = at the disc midplane, and F z — aT 4 , <; = E at the 
photosphere. 



3.2 Grey atmosphere approximation 

For low values of the surface column density, the disc optical 
depth is no longer large. An accurate solution of the vertical 
disc structure would require an excessive amount of com- 
puting time, since frequency-dependent opacities would have 
to be taken into account for consistency. Instead, we chose 
to integrate the set of equations (]l4|-|l5|) down to optically 
thin regions, but we use a grey atmosphere approximation in 
which the temperature varies as T 4 = T s 4 (l/2 + 3/4r), where 
T s is the surface temperature and r the optical depth, given 
by: 

dr 
dz 



Kp. 



(26) 



For optical depths larger than unity, one should use Rosse- 
land opacities kr, whereas at small optical depths, the 
Planck mean opacities ftp are relevant. Here we take 



1 + T. 



1 + Ti 



■Hp, 



(27) 



where r c = Ekr is the estimated disc optical depth. This 
T(t) relation leads to the following boundary condition at 
the surface: 

F z = 2a/" T(t) 4 [E 2 {t) + E 2 (2r s - r)]dr, (28) 
Jo 

where E 2 is the exponential-integral function. This relation 
can be approximated to better than 2% by 
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F z = 2ar s 4 (i 



0.84T, 3/2 e" 2Ta 



(29) 



where r s is the optical depth from the midplane to the sur- 
face. The grey atmosphere approximation is rather crude 
(see e.g. Shaviv & Wehrse 1991; Hubeny 1990); moreover, 
it does not account for energy deposition in optically thin 
layers. However, calculations of the vertical disc structure 
using the Shaviv & Wehrse (1991) radiative transfer code 
show that the grey approximation is quite good , especially 
in vi ew of uncertainties in Q~ mentioned above ([dan et al., 
199S ) . n any case, the grey atmosphere approximation is 
required if one is willing to calculate a large number of such 
disc structures in a reasonable amount of computing time. 

The two photospheric boundary conditions ( poj ) and 
F z = crT 4 are replaced by (E9[) and p vanishingly small (we 
take here p — 1CP 10 g cm~7 The optical depth has to be 
explicitly integrated in the vertical structure. Consequently, 
Eq. ( |2rj ) is added to the set of equations that we solve, with 
the associated boundary condition r = at the disc surface. 



3.3 Numerical method 

We solve the set of equations (|l4|-|l7|) for given values of the 
surface density E and central temperature T c using the same 
adaptive grid numerical method as for the disc equations. 
The grid is now defined by: 



(30) 



+ (e^) + (— ) ( 1+2V >' 

= (^tTc/Q,^) 1 ^ 2 is the vertical scale height of 



where z nolm 
the disc. 

Because we have to solve the vertical structure for any 
given T c and E, a e ff is not known a priori, but is a function 
of T c and E. In other words, there is one more boundary 
condition than the number of differential equations, but one 
parameter is unknown. This can be formulated in a more 
classical problem by introducing a new equation: 



dz 



= 0. 



(31) 



We solve a set of 7 equations |, [T(| |l]) with 7 

boundary conditions in the optically thick approximation, 
while Eq. (^fj|) and the corresponding boundary condition 
are added to the problem in the grey atmosphere approxi- 
mation. We take all weights 0i in Eq. ( |l3| ) equal to 0.5. 



3.4 Results 

The equation of state of matte r is interpolated from the 
tables of Fontaine et al. (1977); in the low temperature 



regime (below 2000 K), which is not covered by these ta- 
bles, Saha equat ions are solved iteratively, as described by 
Paczyhski (1969). The R osseland mean opacities are taken 
from Cox & Tabor (Il976|) above 10,000 K, and from Alexan- 



der (1975)below. The Planck mean opacities are taken from 
Ume ( p.994p , and cover the range 1000 - 30,000 K. The chem- 
ical composition is as sume d to be solar. We note that Liu 
& Meyer-Hofmeister ( 1997 ) have shown that the use of im- 
proved (OPAL) opacities does not much affect the results of 
the vertical structure calculations, since most of the uncer- 
tainty resides in the a-prescription. 

Figure 1 shows an example of the E - T e a "S-curves" 




E (g cm 2 ) 

Figure 1. E - T cff curves for r = 10 10 cm, Mi = 1 M0, and a = 
0.01, 0.1 and 1, in the grey atmosphere (solid line) and optically 
thick (dashed line) approximations. 



we obtained, in both the optically thick and the grey at- 
mosphere approximations. Differences are observed between 
the two cases, essentially when the disc is very optically thin 
(i.e. E less than about 1 g cm -2 ), and when the disc is con- 
vective. The values E max and E m i n which define the upper 
and lower stable branches can be fitted, in the grey atmo- 
sphere approximation, by: 

\ -0.37 

0.85 / Ml 



Emax = 13.1 a 

and 



m u=) ■<»> 



M 



-0.37 



with the corresponding mass transfer rates 

■ \ -0.89 

M" it = 4.0 10 15 a- 0M f Ml 
crlt \^M 

and 



/ r \2.67 



«S,- 8.0 10" a- (^J-)- gs -.(35> 

where M~ it corresponds to E max and M<J it to E m i n . We 
obtain in the optically thick case: 



E max = 13.4 a 



si) iwhrS" «™^ 36 » 



f \ -0.88 2 65 



\ _o ' 89 

M+ 1 = 9.5 10 15 a 01 f g± ) ( 

1 M J VlO 10 cm 



g s 



"(39) 



These are in good agree ment with values obtained previously 
(see e.g. Ludwig, et al. ( 1994 )). A detailed comparison of our 
E — T c ff curves with those obtained using the full radiativ e 
transfer equation is left for a future paper ( [dan et al., 1998 ). 
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Figure 2. Modification of the E — T e g curves when a is no longer 
a constant. The dashed curves correspond to constant a, a = «hot 
(left) and a = a co id (right). The solid line is obtained for a given 
by (fid). The radius of interest is 10 10 cm, and the primary mass 
1.2 M 



4 EVOLUTION OF DWARF NOVA DISCS 

It is well known that models with constant values of a cannot 



reprod uce the observed light curves of dwarf novae (Smak 
1984b) (see however Faulkner et al., 1983). They predict 



that globally, the disc reaches neither the hot nor the cold 
branches. Instead a transition front propagates back and 
forth in a relatively narrow zone, producing rapid and small 
amplitude oscillations of the optical magnitude and the ac- 
cretion rate. One can get around this difficulty if one as- 
sumes that a is different on the hot and cold branch of the 
S-curves. In the following, we assume a = Qhot on the up- 
per branch and a = a co id on the lower branch, with the 
following temperature dependence for a: 



log(a) = log(a co i d ) + [log(a h ot) - log(acoid)] 



1 + 



2.5 x 10 4 K 



.(40) 



The sharp transition between a C oid and evhot is required if 
one wishes to keep the values of E m i n (a) and £ max (cv) in the 
effective S-curve equal to £ m in(ahot) and S max (a co id). 

Figure ^ shows the changes introduced in the E — T e g 
curves when a is given by Eq. (fio]). 

Other a prescriptions, that we do not consider here, 
have been proposed in the literature 



([Meye 



Me 



Hofmcistcr, 1982; |Cannizzo, Chen fc Livio, 1995 ; Vishniac 



& Wheele r 199 6), such as a 



■vc 



ao(H/r) n (see Lasota & 
Hameury ( 1998 ) for a discussion of this prescription). 

In what follows we use a co id ~ 0.01 as it is commonly 
used in the literature. This assumption may be in contra- 
diction with models of tur bulent viscosity in accretion discs 
( |Gammie fc Menou, 1998 ). 



In this article we only consider the case of dwarf no- 
vae in which the disc extends down to the whit e dw arf sur- 
face; its radius is taken from the Nauenberg (1972) mass- 
radius relation for white dwarfs. The visual (V) magnitudes 
are con puted assuming that each annulus of the disk emits 



10 



17 



10 



15 



13 



10 



2.10 



2.05 



2.00 



.o 1.95 
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Figure 3. Outbursts properties for Mi = 0.6 Mq, r in = 8.5 X 10 s 
cm, o cold = 0.04, o hot = 0.20, < r out >= 2 X 10 10 cm, and 
M = 10 16 g s _1 . The upper panel shows the mass accretion rate 
onto the white dwarf, the second one the outer disc radius, the 
third one the disc mass, and the lower panel the visual magnitude. 



blackbody radiation at the corresponding effective temper- 
ature. 

Figures ^| and ^ show two models obtained for Mi = 0.6 



M 



0, n E 



>.5x 10 cm, a co id = 0.04, a hot = 0.20, an average 



and for mass transfer rates of 10 
-i 



g s 



r out of 2 x 10 1U cm, 

(Fig. |) and 10 17 g s _1 (Fig. g). In the first case, the disc 
experiences inside-out outbursts, while in the second case, it 
experiences outside-in outbursts. The results concerning the 
evolution of the disc outer radius are in good agreement, at 
least from a qualitative poin t of v iew, with those previously 
obtained Ichikawa & Osaki (1992). 



The typical surface density and central temperature ra- 
dial profiles that we observe in our simulations are shown in 
Fig. ^. During the inward propagation of the cooling fronts, 
we recognize the density rarefaction wave that precedes the 
region of rapid cooling, while during the propagation of 
the heating fronts, we observe the de nsity spike associated 



with the sharp region of rapid heating ( Pa^ajrjizou^affiknei 



Lin, 1983 



Lin, Papaloizou fc Faulkner, 1985; Cannizzo 
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Figure 4. Same as fig. N. but for M = 10 17 g s 1 . 



1993a). We also note that our light curves are periodic (ex- 
cept when the number of points is so small that numerical 
errors are large). 



4.1 Effect of the outer boundary condition 

We now t ake th e same basic parameters as those used by 
Cannizzo (1993a) in his attempt to reproduce the light curve 
of SS Cyg. We take a h ot = 0.1, a co id = 0.02; the mass 
transfer from the secondary is Mt =10 -9 Mq yr _1 , and 
the mass of the accreting white dwarf is Mi = 1.2 Mq with 
an inner disc radius r; n = Ri — 5 x 10 s cm. We use the 
optically thick approximation for determining the cooling 
function Q~ of the disc. Our reference model uses 800 grid 
points. 

It is worth noting that, because we do not use exactly 
the same cooling function (our values of E m i n and S max 
differ by approximately 10 percent, and we do not use the 
same interpolation procedure between the cold and the hot 
branches of the S-curves, where a varies from a C oid to Qhot), 
we do n ot expect to obtain quite the same result as Cannizzo 
( 1993a ); in particular, our model parameters have not been 
adjusted to fit SS Cyg's observed light curves. 




r (cm) 




r (cm) 

Figure 5. Typical profiles of the surface density S (solid lines) 
and the central temperature T c (dotted lines) observed during the 
evolution of the thin disc, (a) Inward propagation of a cooling 
front and the associated density rarefaction wave, (b) Outward 
propagation of an inside-out heating front and the associated den- 
sity spike. The dashed lines represent £ m i n (upper curve) and 
Smoi (lower curve). 




300 



400 



600 



time (days) 



Figure 6. Calculated light curves assuming that the outer disc 
radius is fixed (lower panel), or that it may vary with time (upper 
panel). The average value of the outer disc radius is the same as 
in the fixed radius case. 



Figure |6| shows the two light curves obtained when the 
outer disc radius is fixed at r out — i x 10 10 cm, and when 
r out is allowed to vary, the constant c in Eq. ^ being such 
that the average value of r out is also 4 x 10 10 cm. 

In the first case, we do reproduce the results of Can- 



nizzo ( 1993a ) qualitatively: we have an alternative sequence 
of short and long outbursts, which are all of the inside-out 
type. There are some quantitative differences, which can be 
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attributed to the differences mentioned above, but the over- 
all aspect is similar. In the second case, we obtain totally 
different results which . not surprisingly, are similar to those 
of Ichikawa & Osaki (1992). We still obtain an alternating 
sequence of long and short outbursts, but both the number 
of short outbursts and the cycle length are quite different. 
As a rule, it is much easier to obtain such a sequence of al- 
ternating short and long outbursts when the outer radius is 
fixed than when it is allowed to vary. 

When the outer edge of the disc is kept fixed at a con- 
stant valu e, we a lways observe i nside -out outbursts, as did 
Canni zzo (1993a), Ludwig et al. (1994) and Ludwig & Meyer 
(1998). This is in contradiction with observations which 
show that, in a number of cases, the instability is triggered 
in the outer parts of the dis c. Suc h a discrepancy has been 
attributed by Ludwig et al. ( 1994] ) to the fact that the disc 
might not extend to the white dwarf surface (note that this 
by itself would not suffice to give outside-in outbursts, since 
the instability would still be triggered in the innermost parts 
of the disc, but it could account for the UV delay). They 
also argued that the viscosity law could be more complex 
that the simple bimodal a prescription. However, when r out 
is allowed to vary, we obtain outside-in outbursts for mass 
transfer rates high enough for the accumulation time of mat- 
ter in the outer disc regions to be shorter than the viscous 
diffusion time. 

The qualitative difference observed between the two 
cases is due to the fact that a much larger fraction of the 
disc is accreted during (large) outbursts when one assumes 
that r out is constant: the density spike associated with the 
heating front brings some material to the outer edge of the 
disc, which splashes on the rigid wall implied by the bound- 
ary condition, and thus strongly increases E; this matter will 
eventually be accreted until the cooling wave starts, when 
E = E m j n at the outer edge of the disc; the accompanying 
rarefaction wave actually brings E below that value. The 
real situation is however quite different: when the density 
spike reaches the disc outer edge, r out increases (matter car- 
rying a large quantity of angular momentum flows outward) , 
and the spike finally dies out. Thus the cooling wave starts 
earlier, at a stage where less mass from the outer disc has 
been accreted; moreover, during the quiescent phase, the 
disc contracts (the outward angular momentum flux is re- 
duced), while mass transfer from the secondary continues. 
Both effects contribute to increase E in the outer disc and 
help the occurrence of outside-in outbursts. 



4.2 Resolution Effects and the Disc Global 
Properties 

The use of adaptive grids is an important improvement on 
previous time-dependent disc models in that the transition 
fronts are always resolved, whatever their location in the 
disc. This is illustrated in Fig. ^where we show the sharpness 
of a heating front and contrast it with the actual width 
"seen" by the adaptive grid. A fair fraction of the grid (~ 100 
of the total 800 grid points in our reference model - adp800) 
is devoted to the front as a region of strong gradients. Note 
that the inner edge of the disc is also resolved by the grid. 

We now evaluate the effect that degrading the numeri- 
cal resolution has on the predicted light curves and on the 
variations of the total disc mass. We compare the relative 




r (cm) 




200 400 

k (grid point) 

Figure 7. An example of the radial profiles of surface density S 
(solid line) and central temperature T c (dotted line) in the real 
space (a) and the grid space (b), using the adaptive grid at a 
resolution N = 800 (reference model). The strong gradients are 
resolved by the grid and about 100 grid points are devoted to the 
resolution of the front itself. 



merits of our code in an adaptive grid version and various 
fixed grid versions that we obtain when the grid function 
W in Eq. ^ is defined as a function of r only. In all cases, 
the outer radius has not been allowed to vary. We tested r, 
log and sqr grids (defined by a linear, a logarithmic and a 
square root spacing between grid points, respectively), but 
we concentrate here only on the results in the sqr grid case 
because this type of grid has been used extensively in the 
literature. Indeed, the viscous equation has a particularly 
simple form when the variab les x — r 1//2 and S = zE are 
used ( |Bath fc Pringle, 198l| ). We report in Table |l| a few 
examples of grids and resolutions that have been used in 
time-dependent studies so far. 

We investigate the fixed sqr and adaptive grids at reso- 
lutions N = 100, 400, 800, 1600 and we refer to these models 
as adplOO to adpl600 and sqrlOO to sqrl600. Figures ^ to 
|ic| show the time evolution of the V magnitude and the to- 
tal disc mass in these models. A straightforward conclusion 
is that a high numerical resolution (N > 400) is required 
to reach a regime where the general shape of the outbursts 
becomes independent of the resolution. This is particularly 
true for the small outbursts, because the general aspect of 
the large ones is obtained at a relatively moderate resolu- 
tion. The adaptive grid naturally requires a smaller number 
of grid points than the fixed grid to reach this regime (typ- 
ically 200 - 400 points, compared with 800 points for the 
fixed grid). It is worth noting that most results in the litera- 
ture have been obtained with resolutions which do not seem 
sufficient to avoid the resolution limits (see Table |l[) . The in- 
adequacy of such a treatment is partially hidden by the fact 
that (1) the viscosity is a free parameter which can be deter- 
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Table 1. A sample of previous time-dependent studies of thin accretion discs, with the type of fixed grid and the numerical resolutions 
used. 'System' refers to a White Dwarf primary (WD, for the dwarf nova studies) and a black hole or neutron star primary (BH, for the 
black hole or neutron star Soft X-ray Transient studies). Note that Ichikawa & Osaki (1992) and Smak (1994b) are the only studies in 
which the outer disc radius is allowed to move. 



Study 



Lin. 



Smak ( ]l984b| ) 
Papaloizou & F aulkn er (1985) 
Mineshige (Il987|) 



Mineshige k, Wheelerj 1289 ) 
Ichikawa & O saki (jl 9?jj| ) 
Cannizzo ( 1993a ") 

Cannizzo ( 1994j ) 

Cannizzo, Chen fc Livi o (1995) 
Cannizzo 



Cannizzo (199 
Ludwig & Meyer (h998|) 




System 


Type of grid 


Num. Resolution 


WD 


sqr 


20-25 


WD 


r/log 


35 


WD 


log 


401 


BH 


log 


21-41 


WD 


r 


35-45 


WD 


sqr 


25-200 


BH/WD 


sqr 


300 


BH 


sqr 


103-1000 


WD 


sqr 


400 


BH 


r/sqr/log 


21-1000 


WD 


? 


200(500) 





Figure 8. V magnitude light curves produced by our numeri- 
cal model, using the adaptive grid at various resolutions: adplOO 
(N=100), adp400 (N=400), adp800 (N=800), adpl600 (N=1600). 
The outer radius is fixed at r ou t = 4 X 10 10 cm. 



Figure 9. V magnitude light curves produced by our numeri- 
cal model, using the fixed sqr grid at various resolutions: sqrlOO 
(N=100), sqr400 (N=400), sqr800 (N=800), sqrl600 (N=1600). 
r ut is not allowed to vary. 



mined only by confronting the time-dependent disc models 
with the observations, and (2) the cooling term Q~ that 
enters the thermal equation is not well known, and varies 
somewhat from one author to another. These two facts may 
account for two differences that we have with Cannizzo's 



(1993a) results, even when we take his assumptions, param- 



eters and number of grid points. Our small outbursts have 
peak amplitudes significantly smaller than the large out- 
burst peak amplitudes. We also find that the evolution of 
the disc mass is almost independent on the resolution of 
the fixed sqr grid, in the sense that ~ 50 % of that mass 
is invariably accreted during large outbursts, whatever the 
resolution. 



4.3 Resolution of the Transition Fronts 

Until recently, numerical simulations were not able to re- 
solve the structure of the transition fronts. For that reason 
it has been very difficult to compare the predictions of semi- 
analytical models with the numerical results, and th e physics 
of the fronts is not yet understood in much detail (Lin , Pa 

1994 Cm 



paloiz- 



Faulkner, 1985; Cannizzo 



199f 



Cannizzo, Chen fc Livio, 19951 |Vishn iac & Wheeler, 199(3 



Vishniac, 1997| |Mcycr, 1984| ; |Mcycr, 1986] ). We leave this 
comparison for a future paper, and here instead we deter- 
mine the numerical resolution required to address this prob- 
lem properly. A heating or cooling front which is not nu- 
merically resolved will be artificially enlarged; the physical 



© 0000 RAS, MNRAS 000, 000-000 



Accretion disc outbursts 11 




aoo 400 
Time (days) 



Figure 10. Evolution of the total disc mass with time for the var- 
ious test models at numerical resolutions TV = 100, 400, 800, 1600. 
The adaptive grid results are shown (solid lines) together with 
the fixed sqr grid results (dashed lines). In all cases, r ou t is fixed. 



results obtained in these conditions are very questionable. 
In the extreme case where a gradient (like a front) is limited 
to one numerical cell (between two grid points), the fluxes 
across the gradient are likely to be reduced and more of a 
numerical than a physical nature. Once again, one should 
keep in mind that uncertainties in the thermal equation (Q), 
and in particular in the radial flux J , may well hide resolu- 
tion effects. For instance, Cannizzo (L993a) shows that the 
light curves of dwarf novae do depend on the precise form 
of the thermal equation. 

We show in Fig. [n] the fractional width, SW/r, of the 
heating and cooling fronts during several outburst cycles 
of our test models. Note that 5W is the width over which 
90 % of the variation of a (defined in Eq. (flc])) occurs. This 
should not to be taken as a definitive value for the front 
width, but rather as a typical width (affected by resolution 
effects as well). What appears clearly in Fig. [H]is that fixed 
sqr grids do not resolve transition fronts well enough in the 
inner parts of the disc, even at high resolution. The discon- 
tinuities seen in the sqr model widths are signatures of the 
resolution limits and are very significant at low resolution 
(TV = 100,400). The situation is systematically worse for 
heating fronts, which are sharper than cooling fronts. The 
adaptive grid has been devised to avoid large temperature 
and density differences between two grid points, and predicts 
the correct width for a moderate number of grid points TV. 

We also note that even when, close to the WD, the 
front widths are limited by the resolution, the resulting light 
curves may still be correct (sqr800-1600). The basic reason 
is that the general outburst shape is not affected by what 
happens close to the compact object, where all characteristic 
timescales are small. On the other hand, an adaptive grid 
may approximately resolve the heating and cooling fronts, 



but produce incorrect results at low resolution (for example 
when TV = 100). This is a result of too many grid points 
being used in steep gradient regions, which leaves too small 
a number of grid points in most parts of the disc. 



4.4 A Comparison of Numerical Grids 

We define the fractional inter-grid spacing (hereafter FIS) 
between two grid points (k and k + 1) as: 

{ Ar \ r k+1 - r k 

V r J k+1/2 



(r fe+1 +r fc )/2' 



(41) 



where the index k refers to the grid point of interest in the 
grid extended from k = 1 at the inner edge of the disc to 
k = TV at the outer edge of the discs. A fixed r grid is defined 
by: 



r k+1 =r k +5, 



(42) 



S being a constant depending on the resolution TV used. In 
this case, the FIS is 



(v 



Ar s 

r / fc+i/2 



ri 1 



N-l n. 
A fixed sqr grid is defined by: 



and the corresponding FIS is 

/ Ar \ _ 2(0^7- y/¥T) 1 
V r ) k+ i/2 ~ N-l 
A log grid is defined by: 



rk+i 



■j3r k . 



The corresponding FIS is 
'A;-\ 2(/?-l) 
1/2 ~ 1 +/3 



(v), 







f r N 
V ri 



(43) 

(44) 

(45) 

(46) 

(47) 
(48) 



The adaptive grid runs show that the heating fronts 
reach a fractional width as low as 8W/r ~ 10~ 2 , and that 
8W/r is ~ 10 -1 ' 3 for the cooling fronts at the inner edge 
of the disc (see Fig. [n]). This means that FISs of the order 



of (HT 1 - 8 , 10-^)/(KT 



10" 



are required if a grid is to 



resolve the cooling and heatin g fr onts with 3/100 grid points 
respectively. We show in Fig. fL2l the FISs for the r, sqr and 
log grids as a function of the numerical resolution TV, and 
we compare them to the FISs required to resolve the cooling 
and heating fronts with 3 and 100 points. The comparison is 
performed at the inner edge of the disc (k — 1) since fronts 
resolved there will be resolved everywhere in the disc: the 
front fractional widths increase with radius while the FISs 
of fixed grids are either constant (log grid) or decrease with 
radius (r and sqr grids). 

The fixed grids require very high numerical resolutions 
to resolve in detail the fronts everywhere in the disc (TV ~ 
10 4 -10 6 ). The situation is worse when the accreting object 
is a black hole or a neutron star: the r and sqr grids will 
resolve the fronts less and less as the inner r adius of the disc 
i s redu ced (see for example Cannizzo et al. (199E), Cannizzo 
The log prescription has the advantage of defining 



a grid with a FIS which does not vary with radius. 
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Figure 11. A comparison of the evolution with time of the cooling and heating front fractional widths (8W/r) as a function of the 
radius of location of the fronts in the disc for several outbursts. The left panel shows models with the adaptive grid at various resolutions, 
while the right panel shows models with the fixed sqr grid at the same resolutions. The inside-out heating fronts are sharper than the 
outside-in cooling fronts, which means that the evolution of the disc is anti-clockwise in these diagrams. The resolution effects are much 
more pronounced in the fixed sqr grid case than in the adaptive grid case and even at high resolution (N = 1600), the fixed sqr grid does 
not resolve heating fronts better than marginally. At lower resolution, the results on the front widths are limited by resolution effects 
and are physically incorrect. At a resolution of N = 100, the front widths are clearly limited by the inter-grid spacing of the fixed sqr 
grid. 



We conclude that reliable studies of the physical proper- 
ties of transition fronts cannot be completed in a reasonable 
amount of computational time by fixed grid codes. An ex- 
ception (see Fig. p"l| ) is the structure of cooling fronts, which 
can probably be addressed, at least in the outer parts of the 
discs, by fixed sqr grid codes at high enough (but reasonable) 
numerical resolutions. 



5 CONCLUSION 

We have constructed a numerical code which can calculate, 
in a reasonable amount of computer time and at very high 
spatial resolutions, long cycles of accretion disc outbursts. 
This code works efficiently in the most general framework 
of the disc instability model and does not require special 
assumptions about viscosity or outer or inner radii. Its va- 
lidity is of course limited by the way physical processes such 
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Figure 12. Dotted-dashed lines show the fractional inter-grid 
spacing (Ar/r, also FIS) required by various fixed grids in order 
to resolve the cooling and heating fronts of our reference model 
with 3 and 100 grid points at the inner edge of the disc (r; n = 5 X 
10 8 cm), and consequently everywhere in the disc (see text). The 
comparison with the theoretical FIS of r (dotted line) , sqr (solid 
line) and log (dashed line) grids, as a function of the numerical 
resolution N, illustrate the very high resolutions needed by fixed 
grids to resolve the transition fronts. 



as turbulent viscosity, convection, radiative transfer etc. are 
treated. Of course it is also a ID code modeling a funda- 
mentally 2D (or even 3D) situation. 

Since the mass of the central object enters the disc equa- 
tions only in Qk , we expect most of our results to be valid for 
BH disc models as well (i.e. similar, but at a slightly smaller 
radius), as long as irradiation and general relativistic effects 
can be neglected. 

In future work we intend to include effects of irradiation 
(Dubus et al. 1998) and to apply the code to a systematic 
study of dwarf nova outbursts and X-ray transient events. 
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